%% Series 3 - Exercise 1 - DOE 2024
%% Experiment on the Young modulus
%
% Table of Contents
% Series 3 - Exercise 1 - DOE 2024
% Experiment on the Young modulus
%     Statement
%     1- Importing data from excel
%     2 - Analyzing of the design
%         Evaluation of the experimental design
%         Replication of the high-low and low-high points
%     3- Analysis of the experiments
%         Calculation of the coefficients of the model with interactions
%             Calculation of the confidence intervals
%             Calculation of the coefficients of the model without interactions
%     4 - Coefficients of the non standardized model
%     5 - Recommendation for future measurements
%%
%
%% Statement
% The data is relative to a set of experiments on steal samples with different concentration in carbon and sulfure, tested at different temperatures in order to determine the elastic behavior of the material.The content of the table corresponds to the experimental conditions as well as the Young modulus following the Hooke phenomenological relation:
% \epsilon=\frac{\sigma}{E}
% with \epsilon the strain, \sigma the longitudinal stress and E the Young modulus.
% Data Loading:
% Load the experimental data from the provided data file into a table for analysis.
% Pre-experiment Model Analysis (before experiment execution):
% Critically analyze the design of the experiment for both a linear model and a model including interactions between carbon (C) and sulfur (S).
% For each model, determine and document:
%             The model equation,
%             The dispersion and correlation matrices,
%             The Variance Inflation Factor (VIF) for multicollinearity,
%             The overall variance of the model.
% Summarize these elements and prepare a set of recommendations for the experimenters, assuming the experiments have not yet been conducted.
% Post-experiment Analysis (after experiment execution):
% Infer the main effects and the interaction between the concentration of carbon (C) and sulfur (S) on the Young’s modulus of the material, using normalized metrics.
% Transform the estimated coefficients back into their original metrics for practical interpretation.
% Confidence Interval Calculation:
% Compute the 95% confidence intervals for all the model coefficients to assess the precision of the estimates.
% Conclusion Preparation:
% Draft a preliminary conclusion for a report that reflects the results of the experiment, focusing on the quality of the measurements and their implications.
% Future Experiment Recommendations:
% Prepare recommendations for a second experimental campaign that includes an evaluation of potential interactions with temperature. Suggest any design improvements needed to study these new factors.
%
%% 1- Importing data from excel

% Setup the Import Options and import the data
opts = spreadsheetImportOptions("NumVariables", 7);

% Specify sheet and range
opts.Sheet = "essais_traction";
opts.DataRange = "A3:G11";

% Specify column names and types
opts.VariableNames = ["Run", "Carbon", "Sulfur", "Temperature", "E1", "E2", "E3"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double"];

% Import the data
data= readtable("DataModuleYoung.xlsx", opts, "UseExcel", false);
disp(data)
% Clear temporary variables
clear opts

%
%% 2 - Analyzing of the design
% The steps are:
% a) Extraction of the matrix of essays from the data

M_essays=data{:,2:4};

% b) Normalization of the matrix of essays

M_essays_st=rescale(M_essays,-1,1,...
    "InputMax",max(M_essays),...
    "InputMin",min(M_essays));
disp(M_essays_st)

% c) Building of the model matrix for the linear model with interaction
%  y=a_o+a_Cx_C+a_Sx_S+a_Tx_T+a_{CS}\ x_Cx_S
% We use the coefficient spec matrix with one column per factor and one row per coefficient.
% Then for our situation the matrix is
% spec = \pmatrix{0 & 0 &0 \cr 1 & 0 & 0 \cr 0 & 1 & 0 \cr 0 & 0 & 1\cr 1 & 1 & 0  }
% for the coefficients a_o, a_1, a_2, a_3 and a_{12}.

spec=[0 0 0;1 0 0;0 1 0;0 0 1;1 1 0];
M_model_st=x2fx(M_essays_st,spec);
disp(M_model_st)

% d) Computation of the dispersion matrix

D=inv(M_model_st'*M_model_st);
list_coef={'a_o','a_C','a_S','a_T','a_{CS}'};
format bank
disp(array2table(D,...
    "RowNames",list_coef,...
    "VariableNames",list_coef));

% You can use the functionsc to visualize the matrix as a color plot:

imagesc(D)
colorbar

% e) Computation of the correlation matrix

C=corrcov(D);
disp(array2table(C,...
    "RowNames",list_coef,...
    "VariableNames",list_coef));

% Wer can observe a correlation of -0.5 between the main effects of carbon and and sulfur, as well as between the constant a_o and the interaction effect a_{SC}.
% A correlation coefficient of -0.5 between the estimates of two main effects in a linear model with interactions suggests a moderate negative correlation. Here are some possible comments and interpretations you could make on this value:
% Moderate Negative Correlation: The value of -0.5 indicates that as one main effect increases, the other tends to decrease, but the relationship is not perfectly linear. It’s a moderate inverse relationship, suggesting that there is some dependency between these two main effects.
% Multicollinearity Concern: While a correlation of -0.5 does not indicate severe multicollinearity (which would generally be observed with correlations closer to -1 or 1), it does point to some level of association between these variables. This could influence the stability of the coefficient estimates in the model. You might want to further explore the Variance Inflation Factor (VIF) for a more robust assessment of multicollinearity.
% Impact on Interpretation: When interpreting the results, it’s important to note that the negative correlation between the two main effects may make it more difficult to separate their individual impacts on the response variable. This interdependence could lead to some confounding in the interpretation of their coefficients.
% In summary, while a correlation of -0.5 is not necessarily alarming, it does warrant attention and further investigation to ensure that the model coefficients are reliable and interpretable.
%
% f) Computation of the variance inflation factors (VIF)

VIF=diag(inv(C));
disp(array2table(VIF,...
    "RowNames",list_coef,...
    "VariableNames",{'VIF'}));

% Theses VIF values between 1 and 1.33 indicate a well-behaved model with respect to multicollinearity. There is minimal correlation between the regressors, and this should allow for clear interpretation of the main effects and interactions. No corrective action is needed for multicollinearity based on these values.
%
%
% g) Plot of the position of the measurement points

fz=14; % fontsize

% plot in 3D
plot3(M_essays_st(:,1),M_essays_st(:,2),M_essays_st(:,3),...
    'ko','MarkerFaceColor','r',"MarkerSize",15)
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)
zlabel('X_T',"FontSize",fz)
box on
grid on
view([-65 40])

% Even if the distribution of the points of measurements is non uniform, the design is spread all over the experimental space. This ensure a good estimation of the linear coefficients (but not for the interaction coefficients between the temperature and the two other factors).
%
%
% h) Computation of the variance function

syms xc xs xt

% model coordinate
fx=[1;xc;xs;xt;xc*xs];

% variance fonction
v(xc,xs,xt)=fx.'*D*fx;

% cursor for changing the third variable
Ref1=0;

% Draw the surfaces for carbon and sulfur
h1=fsurf(xc,xs,subs(v(xc,xs,xt),xt,Ref1),[-1 1]);
h1.MeshDensity=15;
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)
zlabel('var(E)/sigma2',"FontSize",fz)
title('Variance function - linear model with interaction', 'FontSize',14)
% Add the measurement points on the floor 
hold on
plot3(M_essays_st(:,1),M_essays_st(:,2),zeros(9,1),...
    'ko','MarkerFaceColor','r')
hold off
% draw the contour graph
fcontour(subs(v(xc,xs,xt),xt,Ref1),[-1 1],'Fill','on')
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)

% Draw the surfaces for carbon and temperature
h2=fsurf(xc,xt,subs(v(xc,xs,xt),xs,Ref1),[-1 1]);
h2.MeshDensity=15;
xlabel('x_C',"FontSize",fz)
ylabel('x_T',"FontSize",fz)
hold on
plot3(M_essays_st(:,1),M_essays_st(:,3),zeros(9,1),...
    'ko','MarkerFaceColor','r')
hold off

%%
%% Evaluation of the experimental design
% 1. Overall Design Evaluation:
%    The current design is non-optimal, but it remains practicable for the purposes of this experiment. While the design is not ideal, it provides useful data and allows for reasonable conclusions to be drawn.
%
% 2. Dispersion Matrix Analysis:
%    The diagonal elements of the dispersion matrix fall between 14% and 25%, which means that 14% to 25% of the sample variance is being attributed to the estimation of the model coefficients. Ideally, in this case, we would aim for a diagonal element value of 1/9 (approximately 11%). This comparison indicates that the design is somewhat inefficient but still within acceptable limits.
%
% 3. Correlation Between Regressors:
%    The correlation matrix shows that there is moderate collinearity between between the estimators of a_o and a_{CS} and between those of a_Cand a_S. This implies that the estimates of these coefficients are not fully independent, which could slightly affect the precision of the estimated coefficients. However, the collinearity does not reach a critical level that would invalidate the results or require significant corrective measures.
%
% 4. Impact on Confidence Intervals:
%    Due to this moderate collinearity, we expect wider confidence intervals for the affected coefficients. Specifically, the confidence intervals may be inflated by approximately 15% (as the square root of the VIF, 1.33, is 1.15). This increase means that there will be slightly less precision in estimating the effects, but this is still manageable within the experimental constraints.
%
% 5. Variance Function and Sensitivity:
%    The analysis of the variance function suggests that the quality of the model deteriorates when the concentrations of carbon and sulfur move in opposite directions (i.e., one high and one low). This region corresponds to areas where the design lacks coverage, with only a single measurement point in these corners, compared to three measurements in the other corner regions of the design space. This highlights a lack of balance in the design, which affects the model's ability to accurately predict responses in these regions.
%
% 6. Design Improvement Proposal:
%    To address the issue of variance and improve the robustness of the model, it would be beneficial to replicate the experiments at the high-low and low-high points, eventually changing the temperature. Replication in these regions would:
%    - Increase the number of measurements in the areas where the model is currently weakest.
%    - Help reduce variance and improve the stability of the model in those regions.
%    - Likely decrease the confidence intervals by improving the overall balance of the design, leading to more precise estimates of the interaction effects between carbon and sulfur.
%    This adjustment could bring the design closer to the optimal scenario and mitigate the model's sensitivity to specific factor combinations.
%
%% Replication of the high-low and low-high points
% Let's make a proposal to improve the variance function.
% Additional points placed on the plane corresponding to the median temperature ( x_T=0), with the carbon and the sulfure at the oposite extreme high-low and low-high:
%  E_2=\pmatrix{-1 & 0 &0 \cr 1 & -1 & 0 \cr -1 & 0 &0 \cr 1 & -1 & 0  }

New_essais_st=[-1 1 0;1 -1 0;-1 1 0;1 -1 0];

% 1) A new matrix of the model is built, concatenated to the previous one and a new dispersion matrix is computed:

New_model_st=x2fx(New_essais_st,spec); % additional lines of the model matrix
NewD=inv([M_model_st;New_model_st]'*[M_model_st;New_model_st]); % new dispersion matrix

% 2) After that a new variance function V(x_C, x_S, x_T) is built as a symbolic function.

% variance function
New_v(xc,xs,xt)=fx.'*NewD*fx;

% 3) A slider is defined for selecting the value of x_T when representing the variance function as a 2D surface


Ref2=1;

h3=fsurf(xc,xs,subs(New_v(xc,xs,xt),xt,Ref2),[-1 1]);
h3.MeshDensity=15;
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)
zlabel('var(E)/sigma2',"FontSize",fz)

hold on
h4=fsurf(xc,xs,subs(v(xc,xs,xt),xt,Ref2),[-1 1]);
h4.MeshDensity=15;
hold off

% Contour lines of the new variance function are drawn:

fcontour(subs(New_v(xc,xs,xt),xt,Ref2),[-1 1],'Fill','on')
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)

% We observe that replicating the high-low and low-high points would lead to a significant reduction in the variance function at these positions, thereby enhancing the model’s accuracy in these regions. Therefore, this is a strong recommendation to pass on to the experimenters. While this adjustment will not directly improve the model’s resolution with respect to temperature, it will notably enhance the precision of the estimates for the concentrations of Carbon and Sulfur.
%
%
%% 3- Analysis of the experiments
% There are three replicates in the experiment. It is good practice to check whether there are significant differences between these sets. Several methods can be used to assess this:
% Visual inspection: Plot the data in parallel to identify any obvious trends or discrepancies.
% Boxplots: Use boxplots to visually compare the distributions and detect any outliers or shifts in central tendency.
% Statistical comparison: Compare the means and standard deviations of the replicates to quantitatively assess differences and variability.

plot(data{:,5:7},'ro')
xlabel('Runs')
ylabel('Young modulus E')
boxplot(data{:,5:7})
serial_means=mean(data{:,5:7})
serial_std=std(data{:,5:7})

% The plot indicates a similar behavior across the sets. However, the boxplot reveals some minor differences between the three sets. Notably, the second set has a lower average than the others, while the third set exhibits a smaller standard deviation. In this situation, it would be useful to apply a T-test to determine whether these differences are statistically significant:

[h,~,ci]=ttest2(data{:,5},data{:,6})
[h,~,ci]=ttest2(data{:,5},data{:,7})
[h,~,ci]=ttest2(data{:,6},data{:,7})

% Finally, the T-test does not provide sufficient evidence to reject the null hypothesis, indicating that the three sets likely come from populations with the same means and standard deviations. Therefore, we can proceed with confidence, knowing that the three sets can be analyzed together.
%
%
%% Calculation of the coefficients of the model with interactions
% The least squares algorithm will be used to estimate the model parameters. To do this, the model matrix must be constructed to include all three replicates. This is achieved by concatenating the matrix for the basic design three times vertically, ensuring that the replicates are incorporated into the analysis.

% concatenation of the model matrix
Big_model_st=repmat(M_model_st,3,1); 

% formating the vector of response
Response=reshape(data{:,5:7},27,1);         

% least square algoritm
Half_effects_st= (Big_model_st'*Big_model_st)\(Big_model_st'*Response); 

%
% It is interesting to compute the relative half-effects. They ease the analysis letting to seen rapidly the importance of the effects in relation with the constant.
%

Relative_effects=Half_effects_st/Half_effects_st(1);

%
%% Calculation of the confidence intervals

beta=0.95; % the confidence level
alpha=1-beta; % one side risk
t_student=tinv(beta+alpha/2,3*9-5); % student coefficient

residue=Response-Big_model_st*Half_effects_st;
var_exp=residue'*residue/(3*9-5);
D=inv(Big_model_st'*Big_model_st);
CI=t_student*sqrt(diag(D)*var_exp);


Table_effect=array2table([Half_effects_st';CI';Relative_effects'],...
    "RowNames",{'half-effects';' 95% CI';'relative'},...
    "VariableNames",list_coef)
bar(Relative_effects(2:end))
xticklabels({'a_C' 'a_S' 'a_T' 'a_{CS}'})
grid on

hold on
    h5=errorbar(Relative_effects(2:end),CI(2:end)/Half_effects_st(1),...
        "linewidth",6,"LineStyle","none");
hold off

% The effects are small, with each contributing only a few percent. The concentration of carbon has the largest effect, while the interaction between carbon and sulfur is minimal.
% The confidence intervals, which are approximately 1.5, are acceptable for the effect of carbon and reasonably acceptable for sulfur. However, they definitively rule out any significant effect of temperature and the interaction between carbon and sulfur.
% The fit can also be realized with the routine fitlm. Let's do it for the linear model

mdl1=fitlm(repmat(M_essays_st,3,1),Response,"linear")
anova(mdl1,'component',2)
anova(mdl1,'Summary')
plotAdded(mdl1)

%
%% Calculation of the coefficients of the model without interactions
% Let's compute again the coefficient for a linear model without interactions.

M_model_st_lin=x2fx(M_essays_st,"linear");
Big_model_st_lin=repmat(M_model_st_lin,3,1);
Half_effects_st_lin= (Big_model_st_lin'*Big_model_st_lin)\(Big_model_st_lin'*Response);

Relative_effects_lin=Half_effects_st_lin/Half_effects_st_lin(1)

%
%% 4 - Coefficients of the non standardized model
%
% There is two ways to get the non standardized coefficients
% calculate the transformation
% re-compute the effects for the original matrix of essays
% Let's do the transformation with symbolic calculation

syms a_o a_C a_S a_T
syms C S T
syms Cbar Sbar Tbar
syms delta_C delta_S delta_T

x_C=(C-Cbar)/delta_C;
x_S=(S-Sbar)/delta_S;
x_T=(T-Tbar)/delta_T;

y=a_o+a_C*x_C+a_S*x_S+a_T*x_T

y2=collect(y,[C,S,T])

y3=subs(y2,{a_o a_C a_S a_T},{Half_effects_st_lin'});

y4=subs(y3,{Cbar,Sbar,Tbar},{0.05,0.6,0});

y5=subs(y4,{delta_C,delta_S,delta_T},{0.01,0.2,20});

[B,terms]=coeffs(y5)
b=double(B)

%
% Let's proceed with the other method. Let's first compute the model matrix

M_model=x2fx(M_essays,"linear");
Big_model=repmat(M_model,3,1);
Half_effects= (Big_model'*Big_model)\(Big_model'*Response)

% Both methods gives the same results and the model expressed in the units of the laboratory is
% \hat{Y}=184.57+822.47 C -28.71 S -0.06 T
%
%% 5 - Recommendation for future measurements
% The recommendations for future measurements should focus on two key aspects:
% The ability to estimate the desired effects, which depends on the placement of the measurement points.
% The power of the design to detect smaller effects, ensuring sufficient sensitivity.
% In the current design, the measurement points are represented by the red dots. To improve the estimation of interaction coefficients, it is necessary to position the points diagonally across the experimental space, as shown by the green points.

% plot in 3D
plot3(M_essays(:,1),M_essays(:,2),M_essays(:,3),...
    'ko','MarkerFaceColor','r',"MarkerSize",15)
xlabel('x_C',"FontSize",fz)
ylabel('x_S',"FontSize",fz)
zlabel('X_T',"FontSize",fz)
box on
grid on
view([-65 40])
hold on
M_essays_add=[0.04 0.8 20;...
    0.04 0.8 -20;...
    0.06 0.4 20;...
    0.06 0.4 -20;...
    ];

plot3(M_essays_add(:,1),M_essays_add(:,2),M_essays_add(:,3),...
    'ko','MarkerFaceColor','g',"MarkerSize",15)
hold off

% We can check that now the new design allows the estimation of the interaction coefficients:

% standardization
E_new=rescale([M_essays;M_essays_add],-1,1,...
    "InputMax",max(M_essays),...
    "InputMin",min(M_essays));
% the new matrix of the model
M_new=x2fx(E_new,'interaction');
% the new matrix of dispersion
D_new=inv(M_new'*M_new);

disp(D_new)
% the new VIF
VIF_new=array2table(diag(inv(corrcov(D_new))),...
    "RowNames",{'a_o';'a_c';'a_S';'a_T';'a_CS';'a_CT';'a_ST'},...
    "VariableNames",{'VIF'})

% We observe that with the addition of this new point, the VIF reaches 1, which represents the optimal value.
% The next step is to ensure that we can accurately estimate the interaction effect that was previously detected but lacked acceptable confidence. To achieve this, we hypothesize that the sample error will remain at 2.5 kPa, and we aim to determine the number of replicates necessary to make the interaction coefficient statistically significant. Our goal is to target a confidence interval of approximately 0.05 kPa.
% The calculation of the confidence interval is:
% CI_\beta=t_{\alpha/2, \nu} \ \sqrt{ D_{ii} \ s^2}
%

Ncoef=7; % number of coefficients in the model
Nexp=13; % number of experiments per replicates
alpha=0.05; % risk level
Nmax=20; % maximal number of replicates
sample_error=2.5; % sample error
CI_new=zeros(Nmax,1);


for Nrep=1:Nmax
N=Nrep*Nexp; % total number of experiments
nu=N-Ncoef; % degree of freedom
Mrep=repmat(M_new,Nrep,1); % matrix of the model with replicates
Drep=inv(Mrep'*Mrep); % matrix of dispersion with replicates

i=6;
Dii=Drep(i,i); % diagonal element for the coefficient i


CI_new(Nrep)=tinv(1-alpha/2,nu)*sample_error*sqrt(Dii);
end
scatter(1:Nmax,CI_new,200,"blue","filled")
grid on
box on
xlabel('number of replicates')
ylabel('confidence interval (standardized)')

% We observe that with 10 to 12 replicates, it is possible to reduce the confidence interval to the selected target. Therefore, the recommendation is to add the four proposed points to the design and ensure enough replicates are performed to achieve at least 12 degrees of freedom for the residuals.